#!/usr/bin/python
# GalacticModelTemp.py

# Library modules
from time import asctime
from os import linesep
import math
import numpy
#from matplotlib import use as useBackend
#useBackend('Agg')


# Intake parameters
# Units/types noted are the FINAL units/types, after run through alterpars.

# SIMULATION PARAMETERS
import GM_simulationparameters as spars
simpars = {# Parameters directly set by user
           'isacluster' : spars.isacluster,
           'nodes' : spars.nodes,
           'codedir' : spars.codedir,
           'outdir' : spars.outdir,
           'pydir' : spars.pydir,
           #'Nptsxy' : spars.Nptsxy, # int
           'scale' : spars.scale, # pc float
           'galactocentric' : spars.galactocentric, # Not yet implemented
           'LensedMag' : spars.LensedMag, # float # Not yet implemented
           'SourceMag' : spars.SourceMag, # float
           'Magnification' : spars.Magnification, # float # Not yet implemented
           'brightlimit' : spars.brightlimit, # int
           'lres' : spars.lres,
           'bres' : spars.bres,
           'dres' : spars.dres,
           'band' : spars.band,
           'c0_V' : spars.c0_V,
           'rho0_V' : spars.rho0_V,
           'c0_I' : spars.c0_I,
           'rho0_I' : spars.rho0_I,
           'c0' : '',
           'rho0' : '',
           'fitcrho' : spars.fitcrho,
           'fixedrho' : spars.fixedrho,
           'fixedc' : spars.fixedc,
           'latitudeweight' : spars.latitudeweight,
           'latitudeweightexternal' : '',
           'latitudeweightinternal' : '',
           'interpolateMv' : spars.interpolateMv,
           'intervals' : spars.intervals,
           'ApparentLimits' : spars.ApparentLimits,
           'interpolateMass' : spars.interpolateMass,
           'maxmassinterval' : spars.maxmassinterval,
           # Derived parameters
           'ScaleZ' : '', # float pc
           'dlocs' : '', # float index-like array
           'llocs' : '', # float index-like array
           'blocs' : '', # float index-like array
           'delds' : '',
           'dells' : '',
           'delbs' : '',
           'Nptsdlb' : '', # int
           'PtScaled' : '', # float pc
           'PtScalela' : '', # float radians
           'PtScaleba' : '', # float radians
           'timestamp' : '', # str; encodes the time of running
           'codeindex' : -1,
#           '' : spars.,
#           '' : '',
           }

# PHYSICAL PARAMETERS
import GM_physicalparameters as ppars
phypars = {# From Han
           'R0' : ppars.R0, # pc float
           'hR' : ppars.hR, # pc float
           'hdz' : ppars.hdz, # pc float
           'ScaleList' : ppars.ScaleList, # pc float array
           'bulk_vels' : ppars.bulk_vels, #pc/yr float array
           'sigma_vels' : ppars.sigma_vels, #pc/yr float array
           # Derived parameters, Han
           'vanillaScaleList' : '',
           'modifiedScaleList' : '',
           # Binney & Merrifield 1998
           # The luminosity function
           'MvList' : ppars.MvList, # mags array
           'LumFuncList' : ppars.LumFuncList, # stars pc^-3 mags^-1 float array
           # Derived parameters, B & M
           'delMvList' : '',
           'vanillaMvList' : '',
           'vanillaLumFuncList' : '',
           'vanilladelMvList' : '',
           'modifiedMvList' : '',
           'modifiedLumFuncList' : '',
           'modifieddelMvList' : '',
           # Cox 1999
           # The measured angular density of stars
           'bvals' : ppars.bvals, # degrees array
           'MvVals' : ppars.MvVals, # mags array
           'N_bM_V' : ppars.N_bM_V, # float array log10(stars mag**-1 deg**-2)
           'allsky_V' : ppars.allsky_V, # float array log10(stars mag**-1 deg**-2)
           'N_bM_I' : ppars.N_bM_I, # float array log10(stars mag**-1 deg**-2)
           'allsky_I' : ppars.allsky_I, # float array log10(stars mag**-1 deg**-2)
           'N_bM' : '', # float array log10(stars mag**-1 degree**-2)
           'allsky' : '', # float array log10(stars mag**-1 degree**-2)
            # Derived parameters, Cox
           'limitset' : '', # float mags array
           'delMvVals' : '', # float mags array
           'vanillaMvVals' : '', # float mags array
           'vanillalimitset' : '', # float mags array
           'vanilladelMvVals' : '', # float mags array
           'fittingMvVals' : '', # mags array
           'fittinglimitset' : '', # float mags array
           'fittingdelMvVals' : '', # float mags array
           'modifiedMvVals' : '', # float mags array
           'modifiedlimitset' : '', # float mags array
           'modifieddelMvVals' : '', # float mags array
           # Kroupa et al. 1993
           # The mass function
           'masslimits' : ppars.masslimits, # 2d Msun float array
           'massindices' : ppars.massindices, # 1d unitless float array
           'massmults' : ppars.massmults, # 1d float array
           'MvMass' : ppars.MvMass,
           # Derived parameters, Kroupa
           'delMvMass' : '',
           'limitsMvMass' : '',
           'vanillaMvMass' : '', 
           'vanilladelMvMass' : '',
           'vanillalimitsMvMass' : '',
           'modifiedMvMass' : '',
           'modifieddelMvMass' : '',
           'modifiedlimitsMvMass' : '',
           # Derived parameters, Han + Kroupa
           'MvMassScales' : '',
           'modifiedMvMassScales' : '',
           'vanillaMvMassScales' : '',
#           '' : ppars.,
#           '' : '',
           }


# Homemade modules
import GM_alterpars # OK
#import GM_nonlinear
import GM_distinction # OK
import plot_distinction # OK+
import GM_LumFunc # OK
import plot_LumFunc # OK+
import GM_fitc0rho0 # OK
# These modules are all secretly used by GM_fitc0rho0
#import GM_apparent # OK
#import GM_logN # OK
#import GM_N_bvals # OK
# Plot the aftermath, then delete the unloved
import plot_fitc0rho0 # OK+
import plot_apparent # OK+
import plot_logN # OK+
import plot_logN_HanStyle # OK+

import GM_darklight
# Secretly used by GM_darklight and GM_MassDistrib and, well, lots of things
#import GM_InterpRoutines
import plot_darklight
import GM_MassDistrib # OK?
import plot_MassDistrib # OK?+

import GM_velocities # OK?
import plot_velocities # OK?+
import GM_EstimateDependence
import plot_EstimateDependence
import GM_SimpleIntegrate # OK?
import plot_SimpleIntegrate # OK?+
import plot_SimpleIntegrateCluster # OK?+

# Begin
print 'Beginning.'
print asctime()
print linesep

# Update parameters so that their units and datatypes are appropriate
(simpars, phypars) = GM_alterpars.alterpars(simpars, phypars)
del GM_alterpars
#simpars = GM_nonlinear.nonlinear(simpars)
#del GM_nonlinear




#simpars['timestamp'] = '20110430_161047'
#outfilenames = ['/science/afournier/GalacticModel/Output/TempData20110430_161047/Out/Out0.txt', '/science/afournier/GalacticModel/Output/TempData20110430_161047/Out/Out1.txt']
#outfilenames = ['/Users/Amanda/Desktop/TempCopy/Out0.txt','/Users/Amanda/Desktop/TempCopy/Out1.txt']

simpars['timestamp'] = '20110502_145615'
outfilenames = []
name1 = '/Users/Amanda/Desktop/TempCopy/Newest/Out'
name2 = '.txt'
for i in range(0,10):
    name = name1+'0'+str(i)+name2
    outfilenames.append(name)
for i in range(10,68):
    name = name1+str(i)+name2
    outfilenames.append(name)
simpars['outdir'] = '/Users/Amanda/Desktop/TempCopy/Newest/Output/'


### For testing only!
ncodes = len(outfilenames)
### For testing only!
    # Read back in


if ncodes == len(outfilenames):
    import cPickle as cP
#    testfile = open(outfilenames[0])
#    testpars = cP.load(testfile)
#    testfile.close()
#    testApparentLimits = testpars['simpars']['ApparentLimits']
#    testintervals = testpars['simpars']['intervals']
#    test = testpars['simpars']['']
#    test = testpars['simpars']['']
#    import GM_InterpRoutines as IR
#    testMvVals = IR.NewEvenArray(testApparentLimits, testintervals**-1, True, False)
#    print testMvVals
#    phypars['MvVals'] = testMvVals
#    simpars['ApparentLimits'] = testApparentLimits
#    simpars['intervals'] = testintervals
#    simpars[''] = test
    #(lensrates, MvHist, MassHist, DSHist, DLHist, bHist) = ratestuple
#    lensrates
    # Read in and concatenate the results. Mv handled differently here.
    MvHistFull = numpy.zeros(len(phypars['MvVals']))
    MassList = []
    DSList = []
    DLList = []
    bList = []
    for i in range(0, len(outfilenames)):
        outfilei = open(outfilenames[i])
        outpars = cP.load(outfilei)
        outfilei.close()
        theserates = outpars['ratestuple']
        (tlensrates, tMvHist, tMassHist, tDSHist, tDLHist, tbHist) = theserates
        MvHistFull[i] = float(tMvHist[-1]) # Wicked little quirk; quick fix
        MassList.append(tMassHist[-1])
#        print tMassHist
        DSList.append(tDSHist[-1])
        DLList.append(tDLHist[-1])
        bList.append(tbHist[0])
    del theserates, tlensrates, tMvHist, tMassHist, tDSHist, tDLHist, tbHist
    # Arrayify 'em
    MassHistFull = numpy.array(MassList) # 2d
    DSHistFull = numpy.array(DSList) # 2d
    DLHistFull = numpy.array(DLList) # 2d
    bHistFull = numpy.array(bList) # 2d oversize
    # Bring back less-informative versions
    MvValMask = numpy.less_equal(phypars['MvVals'][numpy.newaxis,:],simpars['SourceMag'][:,numpy.newaxis])
    import matplotlib.pyplot as ppt
    def cheat111():
        ppt.subplot(311)
        ppt.subplot(131)
        ppt.subplot(111)
    wallaby = False
    if wallaby == True:
        for i in simpars:
            cmd = "%s = simpars['%s']" % (i,i)
            exec cmd
        for i in phypars:
            cmd = "%s = phypars['%s']" % (i,i)
            exec cmd
        height = math.fabs(MvVals[0] - MvVals[-1])
        cheat111()
        anglelims = [MvMass[0,1], MvMass[-1,1], MvVals[0], MvVals[-1]]
        width = math.fabs(MvMass[0,1] - MvMass[-1,1])
        aratio = width/height
        ppt.imshow(MassHistFull, extent=anglelims, aspect=aratio)
        ppt.title('MassHistFull')
        ppt.savefig(outdir+'MassHistFull.pdf')
        cheat111()
        D = dlocs*PtScaled
        anglelims = [D[0], D[-1], MvVals[0], MvVals[-1]]
        width = math.fabs(D[0] - D[-1])
        aratio = width/height
        ppt.imshow(DSHistFull, extent=anglelims, aspect=aratio)
        ppt.title('DSHistFull')
        ppt.savefig(outdir+'DSHistFull.pdf')
        cheat111()
        anglelims = [D[0], D[-1], MvVals[0], MvVals[-1]]
        width = math.fabs(D[0] - D[-1])
        aratio = width/height
        ppt.imshow(DLHistFull, extent=anglelims, aspect=aratio)
        ppt.title('DLHistFull')
        ppt.savefig(outdir+'DLHistFull.pdf')
        cheat111()
        ba = blocs*PtScaleba
        anglelims = [ba[0], ba[-1], MvVals[0], MvVals[-1]]
        width = math.fabs(ba[0] - ba[-1])
        aratio = width/height
        ppt.imshow(bHistFull, extent=anglelims, aspect=aratio)
        ppt.title('bHistFull')
        ppt.savefig(outdir+'bHistFull.pdf')
#        cheat111()
#        anglelims = [[0], [-1], MvVals[0], MvVals[-1]]
#        width = math.fabs([0] - [-1])
#        aratio = width/height
#        ppt.imshow()
#        ppt.title('')
#        ppt.savefig(outdir+'.pdf')

###
#    MvValMask = MvValMask[:,:ncodes]        
###
    del MassList, DSList, DLList, bList #MvList,
    MvList = []
    MassList = []
    DSList = []
    DLList = []
    bList = []
#    print MvValMask
    for i in range(0, len(simpars['SourceMag'])):
        MaskedMv = numpy.multiply(MvValMask[i][:], MvHistFull)
        MaskedMass = numpy.multiply(MvValMask[i][:,numpy.newaxis], MassHistFull)
        MaskedDS = numpy.multiply(MvValMask[i][:,numpy.newaxis], DSHistFull)
        MaskedDL = numpy.multiply(MvValMask[i][:,numpy.newaxis], DLHistFull)
        Maskedb = numpy.multiply(MvValMask[i][:,numpy.newaxis], bHistFull)
#        Masked = numpy.multiply(MvValMask[i][:,numpy.newaxis], HistFull)
        ReducedMv = numpy.array(MaskedMv)
        ReducedMass = numpy.add.reduce(MaskedMass, axis=0)
        ReducedDS = numpy.add.reduce(MaskedDS, axis=0)
        ReducedDL = numpy.add.reduce(MaskedDL, axis=0)
        Reducedb = numpy.add.reduce(Maskedb, axis=0)
#        Reduced = numpy.add.reduce(Masked, axis=0)
        MvList.append(ReducedMv)
        MassList.append(ReducedMass)
        DSList.append(ReducedDS)
        DLList.append(ReducedDL)
        bList.append(Reducedb)
#        List.append(Reduced)
    MvHist = numpy.array(MvList)
    MassHist = numpy.array(MassList)
    DSHist = numpy.array(DSList)
    DLHist = numpy.array(DLList)
    bHist = numpy.array(bList)
#    Hist = numpy.array(List)
    del MaskedMv, MaskedMass, MaskedDS, MaskedDL, Maskedb
    del ReducedMv, ReducedMass, ReducedDS, ReducedDL, Reducedb
    del MvList, MassList, DSList, DLList, bList
    # Flatten out MvHist to get lensrates
    lensrates = numpy.add.reduce(MvHist, axis=1)
    lensrates = lensrates.flatten()
#    MvHist = numpy.ones((4,len(phypars['MvVals'])))
    ratestuple = (lensrates, MvHist, MassHist, DSHist, DLHist, bHist)
    # Apply a similar process to get the full cumulative map of bHist
    # Outer index is magnitude, from lowest (brightest, least exclusive) to highest (dimmest, most exclusive)
    # so accumulate naturally along this axis
    bHistCumulative = numpy.add.accumulate(bHistFull, axis=0)
    # Inner index is latitude, from lowest (richest, but most exclusive) to highest (poorest, but least exclusive)
    bHistCumulative = numpy.add.accumulate(bHistFull[:,::-1], axis=1)
    CumulativeAngle = simpars['blocs'][:]*simpars['PtScaleba'] # radians
    CumulativeAreaOld = 4*math.pi*(1 - numpy.cos(CumulativeAngle)) # steradians
    CumulativeAngle = CumulativeAngle*(180./math.pi) # degrees
    CumulativeAreaOld = CumulativeAreaOld*(180./math.pi)**2 # sq. degrees
    CumulativeAreaNew = numpy.arange(len(CumulativeAreaOld))
    CumulativeAreaNew = CumulativeAreaNew/float(CumulativeAreaNew[-1])
    CumulativeAreaNew = CumulativeAreaNew*CumulativeAreaOld.max()
    bHistCumulativeAreaOld = numpy.array(bHistCumulative.T)
    import GM_InterpRoutines as IR
    bHistCumulativeAreaNew = IR.OldToNewLinear(CumulativeAreaOld, bHistCumulativeAreaOld, CumulativeAreaNew)
    del IR
    bHistCumulativeAreaNew = numpy.array(bHistCumulativeAreaNew.T)
    bHistTuple = (bHistCumulative, CumulativeAngle, bHistCumulativeAreaNew, CumulativeAreaNew)

    # Pickle all the best results, so that in case something goes wrong with plotting, (or already something went wrong, or I get curious later,) they can be recovered
    resultsDictionary = {'MvHistFull' : MvHistFull,
                         'MassHistFull' : MassHistFull,
                         'DSHistFull' : DSHistFull,
                         'DLHistFull' : DLHistFull,
                         'bHistFull' : bHistFull,
                         'MvHist' : MvHist,
                         'MassHist' : MassHist,
                         'DSHist' : DSHist,
                         'DLHist' : DLHist,
                         'bHist' : bHist,
                         'bHistCumulative' : bHistCumulative,
                         'CumulativeAngle' : CumulativeAngle,
                         'bHistCumulativeAreaOld' : bHistCumulativeAreaOld,
                         'CumulativeAreaOld' : CumulativeAreaOld,
                         'bHistCumulativeAreaNew' : bHistCumulativeAreaNew,
                         'CumulativeAreaNew' : CumulativeAreaNew,
#                         '' : ,
                         }

    outfilename = simpars['outdir']+"FinalTally"+simpars['timestamp']
    outfile = open(outfilename, 'w')
    cP.dump(resultsDictionary, outfile)
    outfile.close()
    del cP



    plot_SimpleIntegrate.plot(simpars, phypars, ratestuple, False)
    plot_SimpleIntegrate.plot(simpars, phypars, ratestuple, True)
    del plot_SimpleIntegrate
    plot_SimpleIntegrateCluster.plot(simpars, phypars, bHistTuple)
    del plot_SimpleIntegrateCluster

#    del bHistCumulative, CumulativeAngle, bHistCumulativeAreaOld, CumulativeAreaOld, bHistCumulativeAreaNew, CumulativeAreaNew
    del MvHist, MassHist, DSHist, DLHist, bHist
    del MvHistFull, MassHistFull, DSHistFull, DLHistFull, bHistFull


#print wiley

print '''That's all for now!'''
print asctime()
print linesep

